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We examine the SPCE and TIP5P water models using a temperature series of MD simulations 
in order to study heat capacity effects associated with the hydrophobic hydration and interaction 
of Xenon particles. The temperature interval between 275 K and 375 K along the 0.1 MPa isobar is 
studied. For all investigated models and state points we calculate the excess chemical potential for 
Xenon employing the Widom particle insertion technique. The solvation enthalpy and excess heat 
capacity is obtained from the temperature dependence of the chemical potentials and, alternatively, 
directly by Ewald summation, as well as a reaction field based method. All three different approaches 
provide consistent results. In addition, the employed reaction field method allows a separation of the 
individual components to the solvation enthalpy into solute/solvent and solvent/solvent parts. We 
find that the solvent /solvent contribution to the excess heat capacity is the dominating contribution, 
being about one order of magnitude larger than the solute/solvent part. The latter contribution 
is found to be due to the enlarged heat capacity of the water molecules in the hydration shell. A 
detailed spacial analysis of the heat capacity of the water molecules around a pair of Xenon particles 
at different separations reveals that the particularly enhanced heat capacity of the water molecules 
in the bisector plane between two Xenon atoms is responsible for the maximum of the heat capac- 
ity observed at the desolvation barrier, recently reported by Shimizu and Chan (J. Am. Chem. 
Soc, 123, 2083-2084 (2001)). The about 60% enlarged heat capacity of water in the concave part 
of the joint Xenon-Xenon hydration shell is found to be the result of a counterplay of strengthened 
hydrogen bonds and an enhanced breaking of hydrogen bonds with increasing temperature. Differ- 
ences between the two models concerning the heat capacity in the Xenon-Xenon contact state are 
attributed to the different water model bulk heat capacities, and to the different spacial extension 
of the structure effect introduced by the hydrophobic particles. Similarities between the different 
states of water in the joint Xenon-Xenon hydration shell and the properties of stretched water are 
discussed. 



I. INTRODUCTION 

Nonpolar simple solutes, such as noble gases or 
methane, don't like to be dissolved in water: They are hy- 
drophobic. Hence the corresponding solvation free energy 
is found to be large and positive [EHil^- However, the 
temperature dependence of the hydration free energy of 
hydrophobic solutes reveals a negative solvation enthalpy, 
which is counterbalanced by an also negative solvation 
entropy P, 0) ■ For the case of the simple solutes the 
latter observation has to be attributed to structural or- 
ganization processes in the solvent phase, due to the lack 
of the solutes internal degrees of freedom. In addition, 
the solvation excess heat capacity is found to be positive, 
typically being in the range between 150 J mol~^ and 
300 J mol~^ around 300 K Q , strongly suggesting an 
enhanced heat capacity of the hydration shell water. This 
essential feature has been used as a key ingredient to 
model hydrophobic hydration effects on the basis of a 
solvent structural (basically "hydrogen bond" ) reorgani- 
zation 4, 5, 6, 7]. 

As a consequence of overlapping hydration shells and 
the release of shell water molecules into the bulk phase. 
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the association of hydrophobic particles is found to be 
strengthened with increasing temperature in order to 
minimize the solvation entropy penalty ^ . This sort 
of entropy-driven association process is usually referred 
to as hydrophobic interaction and has been qualitati vely 
confirmed by numerous simulation studies l a Ifl l hol hU 
111 E III pLIll Oils. 19> 20, 21,|2l|2ll3Ill|2l 
1271 l28l I29I l30| . In line with the negative net entropy 
of association, also a negative heat capacity contribution 
is expected and has been incorporated into models de- 
scribing the hydrophobic association on the basis of the 
solvent exposed surface area ,31,, |3^] . For a full depth re- 
view on the present status of conceptual understanding 
of hydrophobic effects we would like to refer to the re- 
cent review articles of Southall et al. , Pratt [s^l , and 
Widom (33. In addition. Smith and Haymet [s^ provide 
a tutorial overview over related simulation techniques. 

However, recent publications by Shimizu and Chan 
m El Eg, Rick p3 and Paschek Q indicate that 
the heat capacity change, related to the association of 
two hydrophobic particles, exhibits a maximum located 
at distance of the desolvation barrier around 5.5 — 6.5 A. 
The presence and observed strength of the maximum has 
not been considered in models of the hydrophobic inter- 
action based on the solvent accessible surface [HES and 
has also consequences for the heat net capacity change at 
the contact state. Nevertheless, some uncertainty exists 
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about the strength of this effect for hydrophobic parti- 
cles in the contact state. Shimizu et al. find an almost 
vanishing contribution to the heat capacity (nearly 
identical heat capacities for the contact state and the 
fully separated state), whereas Rick 30] observes a neg- 
ative net heat capacity for the associated state, being 
qualitatively in accord with model predictions based on 
the change of the solvent accessible surface. Both results 
were obtained for a pair of methane particles dissolved in 
TIP4P model water. Paschek [s^ finds evidence for both 
scenarios, depending on the water model which is taken 
into account: For the case of the association of Xenon 
particles the SPCE model shows the behavior reported 
by Shimizu et al., whereas the TIP5P model corresponds 
to Rick's scenario. 

In this paper we will provide a quantitative energetical 
and structural study on the origin of the observed heat 
capacity effects. We would like to show that (at least for 
the studied model systems) the solvation heat capacity 
is dominated by the solvent restructuring effects and is 
related to the enhanced breaking of hydrogen bonds in 
different parts of the the hydration shell with increasing 
temperature. In addition, we try to elaborate the rea- 
son for the observed differences between the SPCE and 
TIP5P models. Finally, we would like show that observed 
heat capacity effects can be understood from the proper- 
ties of stretched water. 



the energies obtained from the simulation runs. Since 
this procedure is just intended to show the applicabil- 
ity of our approach, it is applied only to the SPCE 
model. For this purpose we performed a series of sim- 
ulations containing 500 water molecules employing the 
SPCE model ^3 as well as a series containing 500 SPCE 
molecules plus an additional Xenon (Xenon-Xenon pa- 
rameters: a = 3.975 A, efc^^ = 214. 7K 48]) particle. 
The systems were simulated at five different tempera- 
tures 275 K, 300 K, 325 K, 350 K and 375 K at a pressure 
of 0.1 MPa. Each of these simulations extended to 20 ns 
and 2 x 10** configurations were stored for further analy- 
sis. 

To determine the heat capacity change for the asso- 
ciation of Xenon particles we use MD simulations con- 
taining 500 water SPCE and TIP5P .49] molecules and 8 
Xenon particles employing the simulation conditions out- 
lined above. Again, each of the model systems is stud- 
ied at 275 K, 300 K, 325 K, 350 K and 375 K at a pres- 
sure of 0.1 MPa. Here, runs over 60 ns were conducted, 
while storing 7.5 x 10^ configurations for further analy- 
sis. The water/Xenon parameters were obtained apply- 
ing the standard Lorentz-Berthelot mixing rules accord- 
ing to CTij = (cTji -I- ajj) /2 and e^j = ^/eileJJ- 

B. Hydrophobic Hydration and Interaction 



II. METHODS 



MD Simulation details 



We employ molecular dynamics (MD) simulations in 
the NPT ensemble using the Nose-Hoover thermostat 
[37Ll3^ and the Rahman-Parrinello barostat 39,40.] with 
coupling times tt = 1.5 ps and Tp = 2.5 ps (assuming the 
isothermal compressibility to be xt = 4.5 10~^bar~^), 
respectively. The electrostatic interactions are treated 
in the "full potential" approach by the smooth particle 
mesh Ewald summation jlj with a real space cutoff of 
0.9 nm and a mesh spacing of approximately 0.12 nm and 
4th order interpolation. The Ewald convergence factor a 
was set to 3.38 nm~^ (corresponding to a relative accu- 
racy of the Ewald sum of 10~^). A 2.0 fs timestep was 
used for all simulations and the constraints were solved 
using the SETTLE procedure 153]. AH simulations re- 
ported here were carried out using the GROMACS 3.1 
program 13 ■ Statistical errors in the analysis were 
computed using the method of Flyvbjerg and Petersen 
|45| . For all reported systems and different statepoints 
initial equilibration runs of 1 ns length were performed 
using the Berendsen weak coupling scheme for pressure 
and temperature control (tt = Tp = 0.5 ps) [46| . 

In order to obtain the solvation enthalpies and ex- 
cess heat capacities of Xenon, we follow two distinct 
approaches: 1) Indirectly from the temperature depen- 
dence the excess chemical potential. 2) Directly from 



We calculate the excess chemical potential of Xenon in 
water a posteriori from simulation trajectories obtained 
at constant pressure/temperature (NPT- Ensemble) con- 
ditions. For this purpose we employ the Widom particle 
insertion method 50, 51] according to 

^_,^^ {Vfds^+,eM-(3AU) 



IJ-e 



(V) 



where AU = C/(s^+^ L) - U{s^\L) is the potential en- 
ergy of a randomly inserted solute [N -1-1)- particle into 
a configuration containing N solvent molecules. The 
Si — L~^fi (with L = V^l"^ being the length of a hy- 
pothetical cubic box) are the scaled coordinates of the 
particle positions and / sW+i denotes an integration over 
the whole space. The brackets (. . .) indicate isothermal- 
isobaric averaging over the configuration space of the 
A^-particle system (the solvent). The entropic and en- 
thalpic contributions to the excess chemical potential can 
be obtained straightforwardly as temperature derivative 
according to 



dT 



and 



Me 



TSex (2) 



and the isobaric heat capacity contribution according to 



-T 



dT 



(3) 



p \ / p 
In order to perform the calculation most efficiently we 
have made use of the excluded volume map (EVM) tech- 
nique |52l Is^ by mapping the occupied volume onto 
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a grid of approximately 0.2 A mesh-width. Distances 
smaher than 0.7 x cTy with respect to any solute molecule 
(oxygen site) were neglected and the term exp[—[3 AU) 
taken to be zero. With this setup the systematic error 
was estimated to be less than 0.02kJmon^. Although 
the construction of the excluded volume list needs a little 
additional computational effort, this simple scheme im- 
proves the efficiency of the sampling by almost two orders 
of magnitude. For the calculation of the Lennard-Jones 
insertion energies AU we have used cut-off distances of 
9 A in combination with a proper cut-off correction. Each 
configuration has been probed by lO'^ successful inser- 
tions (i.e. insertions into the free volume, contributing 
non- vanishing Boltzmann- factors) . 

We use simulations containing 500 Water molecules 
and 8 Xenon particles to study the temperature depen- 
dence of the association behavior of Xenon. The hy- 
drophobic interaction between the dissolved Xenon parti- 
cles is quantified in terms the profile of free energy (PMF) 
for the association of two Xenon particles. The w{r) is 
obtained by inverting the Xenon-Xenon radial distribu- 
tion functions g(r) according to 



w{r) — —kTlngir) 



(4) 



We use temperature derivatives of quadratic fits of 
w{r, T) to calculate the enthalpic and entropic contribu- 
tions at each Xenon-Xenon separation r. All five temper- 
atures 275 K, 300 K, 325 K, 350 K and 375 K were taken 
into account for the fits . The entropy and enthalpy con- 
tributions are then obtained as 



sir) 



dw{r,T) 
df 



and 



h{r) =^ w{r) + Ts{r) . 



(5) 



(6) 



In addition, the corresponding heat capacity change rel- 
ative to the bulk liquid is available according to 



cp{r) = -T 



d^wjr, T) 
92T 



(7) 



C. "Calorimetric" Analysis 

In order to provide a spacial resolution of the water 
contribution to the solvation excess heat capacity, we 
first calculate the individual potential energies of the wa- 
ter and solute molecules. This is done by a reaction field 
method based on the minimum image convention in com- 
bination with a minimum image "cubic" cutoff. This 
approach has been originally proposed by Neumann [H^l 
and was discussed by Roberts and Schnitker [H^IH^. The 
reaction field approach in general is well suited for our 
purposes since it makes it easy to cleanly according the 
potential energy contributions to individual molecules. 



For convenience we divide the potential energies in con- 
tributions assigned to the individual molecules with 



E 



E, = 




(8) 



where Ei is the potential energy assigned to molecule 
i, M is the total number of molecules. The molecule- 
molecule pair energy 
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(9) 



is then obtained as sum over discrete interaction sites 
Q; and /3, with. Ti^jp — {"^jp ~ ^ia 

I based on the 

molecule/molecule center of mass minimum image sep- 
aration. We employ long range corrections Ei^corr. = 
^tcorr. + ^icorr. accountiug for elcctrostatic, as well as 
Lennard Jones interactions. The electrostatic correction 



E\ 



-d 



27r - 
= — Ddi 
W 



(10) 



is a reaction field term, corresponding to the cubic cutoff, 
assuming an infinitely large dielectric dielectric constant. 
Here rfj = X]q lia^ia is the dipole moment of molecule i, 
D = J2i di is the total dipole moment of all molecules 
in the simulation cell and V is the instantaneous volume 
of the simulation box. Ef^^^^ has also been considered 
as the extrinsic potential and has been shown to provide 
configurational energies quite close to the values obtained 
by Ewald summation (with tin-foil boundary conditions) 
|55| . In order to be consistent with applied cubic cutoff 
procedure for the electrostatic interactions, we also use a 
Lennard Jones correction term for the cubic cutoff 

j p 



E, 



LJ 



V ■ 



63 



h9 {^iajP^ia^ 



(11) 



with 6 = denoting the half box length. Ke = [2 + 

15V2arctan(l/V2)]/6 « 2.5093827 and K12 = [17774-1- 
77409%/2arctan(l/V2)]/207360 « 0.4106497 are analyt- 
ically integrated factors accounting for the cubic cutoff 
geometry. In case of the more common spherical cutoff, 
one would have to replace kq/Ii^ by 47r/3i?3 and Knlh^ 
by 47r/9i?^, with being the cutoff-radius. 

In addition to the procedure outlined section III Bl we 
can as well directly use the individual energies to cal- 
culate the the solvation enthalpies and heat capacities 
according to 

^ex — (-^solute) ^~ 

(j^shou) X ((^shcll) - (£^bulk)) (12) 
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where (-Bsoiute) is the average potential energy of the so- 
lute molecule, (Ssheii) is the average potential energy 
of the water molecules in a sufficiently large solvation 
sphere (here we use a radius of 1.0 nm) around the so- 
lute molecule, whereas (-Ebuik) is the energy of the water 
molecules outside this sphere, (rishcii) is the correspond- 
ing number of water molecules in the solvation sphere. 
From the temperature dependence of /lex we can obtain 
the corresponding heat capacities straightforwardly. 



III. DISCUSSION 
A. Hydrophobic Hydration 

The hydration free energies ^.^x for Xenon, given in Ta- 
ble HI were obtained for a system of 500 SPCE molecules 
employing the Widom particle insertion method. In Fig- 
ure^these data are shown, as well as a quadratic fit of the 
data with respect to the temperature. The corresponding 
solvation entropies and enthalpies are derived from the 
fitted temperature dependence, and are depicted as dot- 
ted and dashed lines in Figure^ respectively. A detailed 
comparison of the thermodynamic solvation properties of 
Xenon and other noble gases using different water mod- 
els, with experimental data, however, has been the sub- 
ject of a previous publication |3^ . In this contribution 
we would like to focus on the energetical and structural 
details regarding the hydration and association of hy- 
drophobic Xenon particles. 

In addition to the indirect determination of the sol- 
vation enthalpy based on the temperature dependence of 
the solvation free energy, also a direct calculation is possi- 
ble. Therefore we have to consider the potential energy of 
the solute molecule, as well as the change of the energy of 
the water molecules in the solvation sphere relative to the 
bulk. Another hrute force approach according to Durell 
and Wallqvist ^5^, however, is to determine the energy 
difference between simulations of two distinct systems, 
one containing only the solvent water, and another, con- 
sisting of water plus an additional solute molecule. This 
procedure has also been applied to the case of Xenon and 



T/K 



^ea;/kJmol 



275 
300 
325 
350 

375 



6.15 ±0.1 
7.71 ±0.1 
9.05 ±0.1 
9.87 ±0.1 
10.41 ± 0.1 



TABLE I: Calculated excess chemical potential for Xenon 
dissolved in SPCE water, obtained from the 500 molecule systems 
by the Widom particle insertion technique. 
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FIG. 1: Squares: Excess chemical potential of Xenon in SPCE wa- 
ter. Solid line: Cubic fit to the excess chemical potential of Xenon 
with respect to the temperature according to fJ,ex{T) = fio + /iiT + 
ti2T'^ with /^o = -35.90 kJmol-\ /ii = 2.33 X lO"! kj mol"! K"! 
and fi2 = -2.94 X 10~*kjmol~^ K~^. Dotted line: Solvation en- 
tropy Sex according to the fitted data. Dashed line: Solvation 
enthalpy hex according to the fitted data. Open circles: Solva- 
tion enthalpies obtained directly by subtracting the total config- 
urational energies from simulations consisting of 500 SPCE -|- 1 
Xenon molecules and simulations containing 500 SPCE molecules. 
The energy data are according to the particle mesh Ewald summa- 
tion. Filled circles: Solvation enthalpies obtained from 500 SPCE 
+ 1 Xenon simulations by using the individually calculated ener- 
gies for the solute, shell water and bulk water. The energies are 
obtained by using the reaction field method discussed in the text. 
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FIG. 2: Solvation enthalpies AE = hex obtained from the 
500 SPCE + 1 Xenon simulations by using the individually 
calculated energies for the solute, shell water and bulk water: 



Solute/Water ' 



and AEwater/Water — "shell X (-Eghell 



^bulk)- The open diamonds denote the Water /Water contri- 
bution, whereas the the open circles indicate the Solute/ Water 
contribution. The black squares indicate the sum of both con- 
tributions. The lines indicate the excess heat capacities: So- 
lute/Water: cp ex = 15.3 J K"-"^ mol"-"^ ; Water/Water: cp ex = 
140.3JK-1 mol-l; Total cp^e^. = 155.6 J K"! mol"! . 



SPCE water. The configurational energies shown in Ta- 
ble^Jwere directly obtained from GROMACS during the 
simulation, employing the particle mesh Ewald summa- 
tion technique (see section lll Al for details concerning the 
setup of the PME). The rather long simulation runs of 
20 ns length provide us sufficiently accurate data to de- 
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(500 SPCE) 


(500 SPCE ± 1 Xe) 




T/K 


{Ec) /kJmor^ 


{Ec) /\i3Tao\-^ 


(ASc) /kJmor^ 


275 


—2407'^ fi + 1 S 4S 1472 + f) 0026"! 


—24086 6 ± 1 


— 13.0 ± 1.6 


300 


-23288.1 ± 0.7 (-46.5762 ± 0.0014) 


-23296.8 ±0.7 


-8.8 ± 1.0 


325 


-22531.2 ± 0.5 (-45.0624 ±0.0010) 


-22536.2 ±0.8 


-5.0 ±0.9 


350 


-21794.3 ± 0.7 (-43.5886 ± 0.0014) 


-21796.0 ±0.6 


-1.7±0.9 


375 


-21071.2 ± 0.8 (-42.1424 ±0.0016) 


-21068.7 ±0.8 


2.5 ± 1.1 



TABLE II: Left columns: Average configurational energies Ec as directly obtained from the simulations of 500 SPCE molecules with and 
without an additional Xenon particle. The values given in brackets denote the energy per molecule. The data corresponds to the Particle 
mesh Ewald summation, applying appropriate cutoff corrections for the Lennard-Jones contributions. Right column: The difference 
between the two energy contributions denotes the solvation enthalpy hex = (AEc) for Xenon in SPCE water. 



T/K 


(-Bsoiutc) /kJmoI ^ 


(-Eshcii) /kJmoI ^ 


(T^-shGll) 


(Sbuik) /kJmor^ 


(AS) /kJ mor^ 


275 


-10.196 ±0.014 


-48.125 ±0.016 


137.50 ±0.04 


-48.102 ±0.007 


-13.5 ±2.4 


300 


-9.846 ±0.014 


-46.527 ±0.014 


135.94 ± 0.03 


-46.527 ±0.006 


-9.9 ±2.0 


325 


-9.453 ±0.013 


-44.980 ±0.013 


133.97 ±0.03 


-45.020 ± 0.005 


-4.1 ± 1.9 


350 


-9.071 ±0.012 


-43.489 ±0.016 


131.28 ±0.04 


-43.545 ± 0.005 


-1.8 ±2.2 


375 


-8.672 ±0.014 


-42.017 ±0.014 


128.18 ±0.03 


-42.099 ± 0.005 


1.9 ± 1.9 



TABLE III: Potential energies of the solute particle and the water molecules obtained by using a cubic minumum image cutoff, a reaction 
field correction and a cubic Lennard-Jones cutoff correction, as outlined in the text. (-Esolutc) is the potential energy of the solute particle, 
(^'shell) is the average potential energy of the water molecules in the spehere of 1.0 nm radius around the solute particle, (ngijon) is the 
average number of water molecule in this sphere. (Sbuik) is the potential energy of the water molecules outside the solute sphere. The 
net solvation enthalpy is obtained as hex = (Ai?) = (Sgoluto) ± ("shell) x ((-E^shell) - (-Bbulk))- 



termine the hydration enthalpies with an errorbar of ap- 
proximately ilkJmol""'^. A comparison of the directly 
obtained h^x data with the hydration enthalpies deter- 
mined from the solvation free-energies is given in Figure 
^ The dashed curve, which is representing the hydration 
enthalpies according to the fitted ^ex{T) data, is mostly 
lying within the errorbars of the directly obtained en- 
thalpies (open symbols). However, the apparent devia- 
tion at 375 K might be attributed to the restriction of 
the fit of the /iea;-data to the finite temperature interval 
between 275 K and 375 K. 

As an alternative to the two simulations approach, we 
calculate the solvation enthalpy based on the "individ- 
ual" potential energies of solute and solvent molecules 
using the reaction field method discussed in section Hi CI 
Table IIIII contains the energies of the solute molecule 
^'solute and the water molecules in the bulk and in the hy- 
dration sphere of radius 1.0 nm, as well as the number of 
the water molecules in the hydration sphere. A compari- 
son of the energy values for the water bulk using the reac- 
tion field with the data for the pure water system accord- 
ing to the Ewald summation (given in Table reveals 
that the reaction field data are lying systematically about 



0.042 kJ moP^ higher than the Ewald data. Although be- 
ing rather small, this noticeable difference is of systematic 
nature and is presumably due to the lack of higher order 
multipole contributions to the reaction field term, which 
are, of course, present in the Ewald calculation [H^. How- 
ever, the differences apparently cancel out when consid- 
ering energy differences between shell and bulk, when 
determining the solvation enthalpies. As shown in Fig- 
ure Q and Tables ITTl and IIIll the obtained data for hf-x is, 
both, consistent with the two simulations solvation en- 
thalpies, as well as the solvation enthalpies obtained from 
the temperature dependence of the solvation free ener- 
gies. Figure[21reports a division of the solvation enthalpy 
hex (total) into contributions according to the solute 
(-^solute) (Solute/Water) and the solvent h^x - (i^soiuto) 
(Water/ Water). The calculated total excess heat ca- 
pacity cp^ex of 155.6 J mol~^ is dominated by the 
Water/Water contribution of 140.3 J K^^ moP^ , whereas 
the Solute/Water part with 15.3 J moP^ contributes 
only to about 10%. The observed value for cp^ex of 
155.6 J K""'^ moP^ , however, does only qualitatively agree 
with the value of about 280 J moP"'^, observed ex- 
perimentally for T = 300K m E3- We would also like 
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FIG. 3: Xenon- Water center of mass pair distribution func- 
tion gxo— w(^) S'iKi relative change of the water potential energy 
Ai?wator(»') = -Bwater(»')-S^^fgj. for all investigated temperatures. 
The arrows indicate the sequence of the g(r)-curves pointing from 
low to high temperatures, a. SPCE model b. TIP5P model. 



to point out that the experimentally observed increase 
of the excess heat capacity with decreasing temperature 
(as discussed in Ref. 01) seems not to be present in 
the SPCE simulation data. However, the change of the 
Xenon solvation enthalpy from negative to positive at 
about 363 K has to be attributed mainly to the potential 
energy change of the of water molecules in the hydropho- 
bic hydration shell. In order to elucidate the spacial (ra- 
dial) extension of the solvent contribution to the excess 
heat capacity, we calculate the potential energy of the 
water molecules as a function of distance to the Xenon 
atom. Figure 01 shows the change of the potential energy 
of the water molecules around a Xenon particle with re- 



spect to the bulk value for the SPCE and TIP5P models 
as a function of temperature. Due to the better statistics, 
the data in Figure |21 were obtained from the simulations 
containing 8 Xenon particles. For completeness, also the 
Xenon- Water center of mass pair distribution functions 
are indicated. A rather strong temperature dependence 
of the potential energy of the water molecules in the first 
hydration shell is clearly evident for both models. A 
rather remarkable observation is, that for the lower tem- 
peratures the water molecules in the distance interval be- 
tween 0.4 nm and 0.5 nm, which are corresponding to the 
first hydration shell, exhibit a potential energy even lower 
than the bulk value. With increasing temperature this 
behavior is reversed and the location of the molecules in 
the hydration shell becomes more and more energetically 
unfavorable. On a qualitative level, both the SPCE and 
TIP5P model, exhibit a similar behavior. The TIP5P 
model, however, shows a more strongly pronounced tem- 
perature dependence, and a more richly structured po- 
tential energy profile. We would like to emphasize that 
for the SPCE model at the lowest temperature we find an 
average potential energy for the water molecules in the 
hydration sphere, which lies below the average bulk value 
(see Table Uni). The fact that DureU and Wallqvist did 
not observe this in their simulations 57] is perhaps due 
to their restricted temperature interval (BOOK to 350 K). 

The temperature dependence of the distance depen- 
dent water potential energies are quantified in terms of 
the configurational heat capacity. Figure 01 shows the 
water heat capacities as a function of the distance to 
the Xenon particle. The data are given relative to the 
value for the water bulk of 60.0 JK""'^ mol~^, obtained 
for the SPCE model, and 90.8JK"^mor\ obtained for 
the TIP5P model. As already evident from Figure O 
the strongest heat capacity effect is observed for the first 
hydration shell. For the TIP5P model a noticeably en- 
hanced heat capacity is also observed for the distance 
range between 0.6 nm and 1.0 nm. Beyond a separation 
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distance of 1.0 nm, however, no significant energy differ- 
ence compared to the bulk can be denoted, indicating a 
properly chosen size of the solvation sphere used for the 
calculation of the hydration enthalpies h^x- 



B. Hydrophobic Interaction 

The hydrophobic interaction between two Xenon par- 
ticles is obtained as a profile of free energy w{r) for the 
association of Xenon particles from the simulations con- 
taining 8 Xenon particles. The w{r), as well as the corre- 
sponding enthalpic and cntropic contributions h{r) and 
—Ts{r) are shown in Figure [S] For both water mod- 
els the "contact state" is defined by the minimum of 
the w{r) function at a distance of about 0.42 nm. With 
increasing temperature the well depth of the minimum 
of the profile of free energy at the contact state drops 
from about -2.21kJmor^ (275 K) to -4.25kJmor^ 
(375 K) for the SPCE model and -1.36kJmor^ (275 K) 
to -4.02kJmor^ (375 K) for the T1P5P model. The ac- 
curacy for each of the w{r) profiles has been determined 
to be about ±0.15 kJmoP^. The contact state minimum 
and the minimum characterizing the solvent separated 
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FIG. 5: Profiles of free energy ui(r) obtained for tlie association 
of two Xenon particles (dashed lines) as well as their enthalpic and 
entropic contributions w{r) = h{r) — Ts{r) for all five temperatures 
(275 K, 300 K, 325 K, 350 K, 375 K). Solid lines: h{r). Dotted 
lines; -Ts{r) Top: SPCE model; Bottom: TIP5P model. The 
arrows indicate the sequence of the data curves pointing from low 
to high temperatures. The data shown here is identical to the data 
presented in Ref. I36l . 
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FIG. 6: Relative change in heat capcacity Acp{r) = cp{r) — c^p^^ 
for the hydrophobic interaction between two dissolved Xenon par- 
ticles using the SPCE and TIP5P models. The data correspond to 
the temperature derivative of quadratic fits of ui(r,T) obtained for 
300 K. The solid lines are just to guide the eye. The data shown 
here is identical to the data reported in Ref. 
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FIG. 7: Schematic illustration of the employed definition for the 
cylindrical distribution of the water density and heat capacity 
around a pair of Xenon particles obtained for different Xenon- 
Xenon separation distances dxe-Xe- 



state, which is located at a distance of about 0.72 nm 
to 0.78 nm, are separated by the so called desolvation 
barrier, which is found approximately at a distance of 
0.6 nm. Figure indicates that for lower temperatures 
the contact state is entropically stabilized and enthalpi- 
cally destabilized, in accordance with the observation of 
Smith and Haymet [ll^ and others HI IH HI IH ^ . 
For the TIP5P model, we observe a strong temperature 
dependence of the curves, showing a decrease of the ab- 
solute values for the enthalpy and entropy contributions. 
In contrast to the TIP5P model, the SPCE model does 
exhibit only a weak temperature dependence of h{r) and 
—Ts{r) at the contact state . At the desolvation barrier, 
however, the inverse behavior is observed for both mod- 
els: with increasing temperature the entropy/enthalpy 
compensation effect is found to be enlarged. 

The observed temperature dependence of the enthalpy 
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FIG. 8: Cylindrical distribution function of the SPCE water molecules center of mass around a pair of Xenon particles gxc— Xc— w{^"i 
obtained for a certain Xe-Xe distance interval at T = 300 K. Left: "contact" (dxc— Xc < 0.45 nm). Middle: "desolvation barrier" (0.55 nm< 
dxe— Xe ^ 0.65 nm). Right: "separated" (0.9 nm < dxe— Xe 1-0 nm). 
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FIG. 9: Cylindrical distribution function of the configurational contribution to the heat capacity of the SPCE water molecules around a 
pair of Xenon particles cp(r, z) obtained for a certain Xe-Xe distance interval. The cp{r, z) data have been obtained as linear fits to the 
corresponding water potential energies i?water (''i ^) for all five temperatures. The average bulk value is 62.1 J K~^mol~^ . Left: "contact" 
(dxe— Xe < 0-45 nm). Middle: "desolvation barrier" (0.55 nm < dxe—Xe < 0-65 nm). Right: "separated" (0.9 nm < dxe— Xe 1-0 nm). 



profiles is quantified in terms of a heat capacity profile for 
the association of two Xenon particles, shown in Figure 
|S| In accordance with the observations of Shimizu and 
Chan [23, 115 , and Rick ^3 j we observe a maximum 
of the heat capacity located at the desolvation barrier. 
For the contact state, however, the two water models 
employed in our study show a quite different behavior. 



For the SPCE model we find almost no difference in the 
heat capacity between the contact state and the com- 
pletely separated state, corresponding to the observation 
of Shimizu and Chan jSlll for Methane in TIP4P wa- 
ter. For the TIP5P model, however, we observe a consid- 
erable negative heat capacity contribution at the contact 
state. 
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FIG. 10: Distribution function of the SPCE water molecules center 
of mass around a pair of Xenon particles in the Xe-Xe bisector plane 
gxe— Xe— w(^i ^ = 0) obtained for a certain Xe-Xe distance interval 
at T = 300K. "contact": dxe— Xe < 0.45 nm. "desolvation barrier": 
0.55 nm < dxo— Xe ^ 0.65 nm. "separated": 0.9 nm < dxo— Xc ^ 
1.0 nm. For comparison, also the Xe-W pair distribution functions 
are given. Top: SPCE model. Bottom: TIP5P model. 
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FIG. 11: Heat capacity of the SPCE water molecules around a 
pair of Xenon particles in the Xe-Xe bisector plane cp (r, z = 0) 
obtained for a certain Xe-Xe distance interval. The cp(r, z = 0) 
data have been obtained as linear fits to the corresponding water 
potential energies Ewate.r(T, z = 0) using all five temperatures, 
"contact": dxo— Xo < 0.45 nm. "desolvation barrier": 0.55 nm < 
<^Xc-Xc 0.65 nm. "separated": 0.9 nm < dxe— Xc 1-0 nm. For 
comparison, also the Acp(r) + cp buii; data of Figure 111 are given. 
Top: SPCE model. Bottom: TIP5P model. 



At a first glance this behavior can be rationahzed as 
a simple consequence of the overlapping of the hydra- 
tion shells in case two hydrophobic particles approach 
each other. The positive association enthalpy and en- 
tropy corresponds to a picture of the release of shell water 
molecules with a lower potential energy and higher order- 
ing into the bulk phase. However, when taking the tem- 
perature dependence of the h{r) and s(r) profiles into ac- 
count, two observations are not consistent with this sim- 
ple model considerations: First, since the hydration shell 
water molecules have been shown to exhibit an enhanced 
heat capacity, a reduced number of water molecules in 
the joint hydration shell of a pair of Xenon atoms in the 
contact state should also lead to a net negative heat ca- 
pacity. This, however, is apparently not the case for the 
SPCE model. The second obvious inconsistency is, of 
course, the enlarged positive heat capacity observed at 
desolvation barrier, found for both water models here, 
and even more water models, as discussed in Ref. [s^ . 

In order to elucidate the observed heat capacity pro- 
file for the Xenon- Xenon association, we calculate two- 
dimensional distribution functions describing the water 
density as well as the water heat capacity around a pair 
of Xenon particles, found at a certain distance dxc-Xo- 



These distributions are shown in Figures |S1 and |51 Figure 
13 illustrates the definition of distribution functions. We 
would like to point out that a conceptually similar plot 
of the molecule densities around a pair of hydrophobic 
Methaiie particles has been recently reported by Gosh et 
al. [2^123 ■ If a pair of Xenon atoms is found to belong to 
a certain distance dxe-Xe interval corresponding to one 
of the states "contact" , "desolvation barrier" , or "sepa- 
rated" (see the figure captions of Figures |H1 and El for the 
corresponding definition intervals), the properties of the 
surrounding water molecules (normalized density, poten- 
tial energy) are mapped with respect to their cylinder 
coordinates r and z. z = is chosen in such a way that 
it indicates the bisector plane between the two adjacent 
Xenon particles. 

Figure IHl shows the normalized cylindrical distribution 
functions of the water molecule density around a pair of 
Xenon atoms (7xo-Xo-w(^i 2)1 obtained from the SPCE- 
Xenon simulations at 300 K. In order to provide a more 
quantitative comparison of the TIP5P and the SPCE 
models, we show in Figure [TUl the distribution functions 
obtained for the bisector plane 5xe-Xe-w(''; z = 0). Fig- 
ure |H1 reveals an enhanced water molecule density in the 
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range where the hydration shells of the two Xenon par- 
ticles overlap. With dxe-Xc the gxe-Xo-w(^i z — 0) 
function approaches the 5xc-w('') pair correlation func- 
tion. With increasing dxc-Xc the first peak starts to in- 
crease showing roughly a doubled height (with respect to 
the first peak of gxc-wir)) of the first maximum at the 
desolvation barrier. A further increase of dxc-Xc leads 
again to a decrease in height of the first peak. This is 
qualitatively in accord with the observations of Ghosh 
et al. [13 obtained for Methane in TIP3P water. We 
find that both water model models show qualitatively 
the same behavior. The TIP5P model, however, exhibits 
a slightly stronger increase of the first peak at the desol- 
vation barrier. 

The heat capacity of the SPCE water around a pair of 
Xenon particles cp(r,z) is shown in Figure El The data 
were obtained by sampling the water potential energies 
as a function of r and z for all temperatures and sub- 
sequent linear fitting of the average values with respect 
to the temperature. The potential energies of the wa- 
ter molecules were obtained the same way as discussed 
in the previous sections. The color coding in Figure |2| 
is chosen in such a way that the lower limit coincides 
with the average bulk value for waters configurational 
heat capacity of 62.1 JK""'^ mol"""^. Since we consider here 
a system consisting of 8 Xenon particles and 500 water 
molecules, the bulk value for cp of water is slightly larger 
than the value obtained for pure water discussed in sec- 
tion llll A1 The color spectrum represents a total range of 
30JK"^mor\ Fi gure 121 indicates that virtually no wa- 
ter molecules are found at distances smaller than about 
0.3 nm to any Xenon site, hence no water cp data is avail- 
able. Therefore this region is represented in Figure 
by the color indicating the lower limit for Cp. In order 
to provide a more quantitative representation. Figure ITTl 
shows the heat capacity of the SPCE and TIP5P water 
molecules in the bisector plane. 

Figure El reveals the origin the heat capacity maximum 
at the desolvation barrier, shown in Figure El The wa- 
ter molecules adsorbed close to the bisector plane be- 
tween the two hydrophobic particles exhibit an about 
60 % increased heat capacity (in the "desolvation bar- 
rier" state) , which apparently overcompensates the effect 
of the shrinking of the total hydration shell. In addition, 
Figure El shows that the increased heat capacity of the 
hydration shell water molecules located in the bisector 
plane persists in the contact state, still partially com- 
pensating the reduced solvent accessible surface. The 
differences observed for the contact state of the TIP5P 
and SPCE models seem to be related to the individual 
strengths of these compensating effects. Figure indi- 
cates that the SPCE and TIP5P models behave qualita- 
tively similar. Even the relative changes are of the same 
size: The maximum at the desolvation barrier shows an 
increase of the heat capacity of about 40 J mol~^ for 
both water models. However, the TIP5P water exhibits 
an about 30 JK~^ mol~^ larger bulk heat capacity. As a 
consequence, the effect due to the decrease of the solvent 



accessible surface has to be stronger in case of TIP5P, 
which might explain the negative heat capacity at the 
contact state. In addition, the heat capacity effect might 
as well be influenced by the second hydration shell of 
Xenon (see Figure . The heat capacity of the water 
in the second hydration shell is clearly more strongly af- 
fected in case of the TIP5P model. The contribution of 
this extended hydration shell should also lead to a more 
negative net excess heat capacity. The observed water 
heat capacities provide an explanation for the maximum 
of the heat capacity profile of two associating hydropho- 
bic particles. Hence the differences between the SPCE 
and TIP5P models seem to be related to the delicate 
balance of two compensating effects: the reduction of 
the solvent accessible surface and the increased heat ca- 
pacity of the water molecules in the concave part of the 
joint Xenon-Xenon hydration shell. 



C. Hydrogen Bonding 

Finally we would like to discuss the origin of the ob- 
served heat capacity effects. Therefore we calculate the 
(binding) molecule pair energies impair between water 
molecules which are potentially hydrogen bonded in the 
sense that their 0-0 distance is smaller than 0.35 nm 
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E . / kJ mol 
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FIG. 12: Bottom: Probability density of observing a water- 
water pair energy of water molecules which are separated less than 
0.35 nm, and which are located in the "BULK" phase. Top: Dif- 
ference between the corresponding probabiliy densities of water 
molecules belonging to the Xenon hydration shell ("SHELL"), as 
well as close to the bisector plane between two Xenon particles lo- 
cated at the desolvation barrier ("DSBAR"), and the probability 
density obtained for the bulk (exact definitions are given in the 
text). The arrows indicate the sequence of curves pointing from 
low to high temperatures. 
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FIG. 13: Top: Average number of nerarest (water) neighbors of 
of a SPCE molecule. Bottom: Average nearest neighbor pair inter- 
action energies of a SPCE water molecule belonging to the states 
"BULK", "SHELL", and "DSBAR" (Definitions are given in the 
text). 
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FIG. 14: Potential energy of a water molecule according to near- 
est neighbor interactions (applying an oxygen-oxygen cutoff of 
rc = 0.35nm) belonging to the states "BULK", "SHELL", and "DS- 
BAR". 



|lCt l61| . For this purpose we distinguish between three 
distinct states: "BULK": The separation from the hy- 
drophobic particle is larger than 1.0 nm. "SHELL" : The 
water molecules are found to lie within a sphere with 
radius 0.55 nm around the hydrophobic particle. "DS- 
BAR" : The water molecules are found in the region close 
to the bisector plane between two adjacent Xenon parti- 
cles found within the distance interval between 0.55 nm< 
dxe-Xc < 0.65 nm, in the region defined by r < 0.4 nm and 
— 0.15nm < z < 0.15nm (see Figures [7| and |H1 for further 
explanations). For the first two categories we employ 
simulations containing a single Xenon atom, whereas for 
the "DSBAR" category we use the simulations with 8 
Xenon particles. We would like to point out that the 



simulations containing 8 Xenon atoms provide qualita- 
tively similar results for the first two categories. All data 
reported in this section were obtained using the SPCE 
model, but the TIP5P model is found to furnish a quali- 
tatively similar picture. 

In Figure 1121 the probability densities for obtaining a 
certain water-water pair energy P(£'pair) are given. Com- 
paring the "SHELL" state with the "BULK" state we 
find an increased population of states with binding en- 
ergies less then approximately — 12kJmon^. Interest- 
ingly, the water molecules in the "DSBAR" state show 
an even more enhanced population of strongly bound 
water molecules. A noteworthy observation is that the 
temperature variation of the probability density distribu- 
tions of the "SHELL" and "DSBAR" states with respect 
to the bulk is found to be only small. At least, signifi- 
cantly smaller than bringing a water molecule in either 
the "SHELL" or "DSBAR" state. Figure US contains 
the obtained average pair energies, indicating that the 
binding energies according to the "SHELL" and "DS- 
BAR" states are shifted to lower values. A consequence 
of the small temperature variation of the water-water 
probability density is that the energy averages (impair) 
for the three different states are running mostly par- 
allel to the bulk data. A second parameter control- 
ling the potential energy of the water molecules is, of 
course, the number of possible binding partners, or num- 
ber of nearest neighbors (A^ncighbor) • Figure El there- 
fore shows the number of nearest neighbors as a func- 
tion of temperature. To obtain a first order approxi- 
mation to the water potential energy we may just con- 
sider the first solvation shell of each water molecule, 
according to £; « l/2{Epi,ir) x (iV„eighbor) • The cor- 
responding energy data are shown in Figure 1141 The 
temperature dependence reveals that the water heat ca- 
pacity effects discussed in sections IIII Al and IIII Bl have 
to be largely attributed to the waters nearest environ- 
ment. The heat capacities according to the data shown 
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in Figure d are 52.1JK-^mor\ 58.9 J K"^ mor\ 
and 75.5JK-Vor^ for "BULK", "SHELL", and "DS- 
BAR" , respectively, hence accounting for about 85 % of 
the heat capacity considering the fuU environment. We 
would like to stress two observations: The reduction of 
the number of water neighbors, i.e. reduction of the lo- 
cal water density, leads apparently to a strengthening 
of the water/water pair interactions. We would like to 
emphasize that this behavior is analogous to the behav- 
ior observed for stretched water |3 |6^ , as 
well as to the behavior found in low density patches due 
to density fluctuations of water at ambient conditions 
mill. In Ref. [13 and ^ Sciortino at al. showed 
that a decrease in density is accompanied by a strength- 
ened water-water binding energy. As mechanism they 
identify a decreasing amount of "fifth neighbor" config- 
urations, contributing energetically unfavorable "bifur- 
cated" hydrogen bonds. The same behavior is observed 
in the hydrophobic hydration shell. Hence the water in 
the "SHELL" and "DSBAR" states might be considered 
as locally stretched water, while representing different de- 
grees of stretching. The increased heat capacity with 
stretching is also consistent with the experimentally ob- 
tained heat capacity of pure water based on the lAPWS- 
95 formulation according to Wagner and Prufi jH3|.The 
density dependence of the heat capacities are shown in 
Figure [T51 for several temperatures. As a rough estimate 
we obtain {dcp{p)/dp)T — —5.5 x 10^ kJ m^ kg^^ K^^ for 
T = 300K around p— 1000 kg m"'^. The configurational 
contribution to the heat capacity might by approximated 
as cp{liq.) — cp{gas) ~ 2.27kJkg~^ K~^. Hence a den- 
sity decrease of about 14%, as found for the "DSBAR" 
state, should lead to an increase of the configurational 
heat capacity of about 34%. A density decrease of about 
7%, as observed for the "SHELL" region, should result 
in a heat capacity increase of about 17%. Both are quite 
close to the observed simulation data. In accordance with 
recent ideas of Ashbaugh et al. j63|, the effects of the hy- 
drophobic hydration and interaction discussed here might 
therefore simply refiect features of waters unique equa- 
tion of state. 

Our simulations indicate that the water-water bind- 
ing energy-decrease counterbalances (and at lower tem- 
perature even overcompensates) the effect of the dilu- 
tion of water molecules in the vicinity of hydrophobic 
particles. With increasing temperature the apparent en- 
hanced disintegration of the water network (number of 
nearest neighbors in FigureEJ close to hydrophobic par- 
ticles is therefore responsible for the observed heat capac- 
ity effects. To illustrate this, we finally discuss hydro- 
gen bonding based on a topological criterion. A pair of 
water molecules is considered to by "hydrogen bonded" 
when the 0-0 distance is smaller than 0.35 nm and the 
H - O • • • O angle is smaller than 30° jH [^. Each wa- 
ter hydrogen is considered to be a hydrogen bond donor 
and Jhb denotes the fraction of hydrogen bond donors 
which are not involved in a hydrogen bond, or the frac- 
tion of broken hydrogen bonds, /hb as a function of 
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FIG. 16: Fraction of broken hydrogen bonds ^hb for water 
molecules belonging to the states "BULK", "SHELL", and "DS- 
BAR". 

temperature is shown in Figure 1161 for water molecules 
belonging to different states. In parallel with observed 
increasing potential energies, the increasing temperature 
leads to an enhanced breaking of hydrogen bonds. This 
effect is found to be enhanced in the hydration shell of the 
hydrophobic particles and particularly strong for water 
molecules in the "DSBAR" state. 



IV. CONCLUSIONS 

We have shown that for simple hydrophobic solutes 
the observed positive excess heat capacity of solvation is 
largely (about 90%) determined by an increase of the heat 
capacity of the water in the first hydration shell. The ef- 
fect is mainly attributed to the altered potential energy 
of the water molecules in the first hydration shell, chang- 
ing their energy state from "attractive" at lower temper- 
atures to "repulsive" at higher temperatures. The two 
model systems studied, SPCE and TIP5P water, show 
a qualitatively similar behavior, although the effect is 
found to be quantitatively stronger in case of the TIP5P 
model. 

Our model calculations, in accordance with the calcu- 
lations of Shimizu and Chan [2^ , as well as Rick 
and Paschek [sB] i indicate that the heat capacity for 
the association of two hydrophobic particles exhibits a 
maximum located at the desolvation barrier. Our sim- 
ulations indicate that this seems to be a characteristic 
feature of the hydrophobic interaction of small apolar 
particles. The observed behavior is found to be a coun- 
terplay of two effects: The reduced solvent accessible sur- 
face, and the increased heat capacity of water molecules 
located in the gap between the two hydrophobic parti- 
cles. In case of the SPCE model, the water molecules 
located in this bridging situation show a heat capacity 
increase of about 60 % compared to the bulk, apparently 
overcompensating the effect of a reduced solvent acces- 
sible surface. The differences found for the SPCE and 
TIP5P models in the contact state are accordingly re- 
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lated to quantitative differences in strength of this two 
compensating effects. 

A detailed analysis of the water- water pair interactions 
in the different states (bulk, hydration shell, joint hy- 
dration shell of two Xenon atoms) reveals that the ob- 
served heat capacity effects can be rationalized as a coun- 
terbalance of strengthened hydrogen bonds and an en- 
hanced disintegration of the hydrogen bond network with 
increasing temperature. The reduced number of water 
neighbors in the different parts of the joint hydrophobic 
hydration shell might be interpreted as a locally stretched 
water/hydrogen bond network. A comparison with ex- 
perimental data indicates that the observed heat capacity 
effects have the same order of magnitude as it would be 



expected for pure water with an approximately equally 
reduced density. 
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